% -----------------------------------------------------------------
%   1-D, 1-G DIFFUSION-BASED, NONLINEAR RESPONSE MATRIX SOLVER
% -----------------------------------------------------------------
%  j. roberts, 12/29/2009

% to do: 
%   -- make geometry easier to deal with
%   -- put newton solver in own function
%   -- plotting stuff in own function

clear; clc;
format long

%profile on

% declare global problem data
global ne dc sa ns sw bcL bcR gm pre

problem = 1; ne=100;
seed    = 1;  % seed the nonlinear solvers with one r-b iteration
pre     = 2;  % 0 = none, 1 = jacobi, 2 = appx. jacobian

%----------------------- GET PROBLEM DATA -------------------------
[ne,dc,sa,ns,sw,bcL,bcR,k,j] = probdata(problem,ne);
j       = j/sum(j); % normalize currents & form vector of unknowns
x       = [j k]';
tol     = [1.e-11 1.e-11]; % tolerance for norms


%---------------------- RED BLACK SOLVER --------------------------
disp('----RED BLACK----')
mxrbit = 20;
tic
[x1,itrb,nrmrb,xx1] = redblack(ne,dc,sa,ns,sw,x,tol,mxrbit);
t1 = toc;
plotrb = nrmrb/nrmrb(1);  % relative nonlinear residual, for plot
xsol = xx1(:,end);        % the "solution"
p1 = convest(xx1,xsol,1); % estimate rate of convergence
n1=norm( feval('respfct',x1) );  % compute the final norm
xsol    = [xx1(:,end)' 1.0]';    % save the benchmark solution

%profile on
%---------------------- NEWTON-GMRES SOLVER -----------------------
% uses newton-gmres functions written by kelley
disp('----NEWTON-GMRES----')
tic
if seed==1
    [x,tmp,tmp,tmp] = redblack(ne,dc,sa,ns,sw,x,tol,1);
    j = x(1:end-1);
    x(1:end-1) = j/sqrt(j'*j); % pre-normalize 
end
xo = [x' 1.0]'; 
nrm0 = 1;
[x2, ithistGM, ierr, nrmgm, xx2] = nsolgm2(xo,'respfctAUG',tol,nrm0);
if nrmgm(end) == 0
    nrmgm(end) = eps; % in case the last delta is ~ 0.
end
t2 = toc;
itgm = length(ithistGM);
disp(['GMRES per NEWTON = ',num2str( sum(ithistGM(:,3)/(itgm-1)))])
plotgm = ithistGM(:,1)/ithistGM(1,1); % relative nonlinear residual
p2 = convest(xx2,xsol,0); % est. rate of conv. w/ RB soln. as ref
%profile viewer
%---------------------- NEWTON SOLVER -----------------------------
% uses the standard newton method w/ a fd-approximated jacobian
disp('----NEWTON----')
gm = 0;
tic
if seed==1
    [x,tmp,tmp,tmp] = redblack(ne,dc,sa,ns,sw,x,tol,1);
    x(1:end-1) = j/sqrt(j'*j);
end
x3          = [x' 1.0]'; 
xx3(:,1)    = x3;
z           = respfctAUG(x3); % initial residual
itmx        = 25; 
it          = 1;
zz(1)       = norm(z);
nrmcn(1)    = 1;
while ( zz(it) > tol(1)*zz(1)+tol(2) && it < itmx )
    fp          = jacob('respfctAUG',x3,z);
    s           = -fp\z;
    x3          = x3 + s;
    it          = it+1; 
    xx3(:,it)   = x3; % Keep all x's to estimate rho
    z           = respfctAUG(x3);
    zz(it)      = norm(z);  
end
t3 = toc;
plotcn = zz/zz(1); % relative nonlinear residual
p3 = convest(xx3,xsol,0); % est. rate of conv. w/ RB soln. as ref

%--- OUTPUT
semilogy(0:itrb-1,plotrb,'k',0:itgm-1,plotgm,'b--',...
         0:it-1,plotcn,'r-.')
axis([0,max(itgm,max(itrb,it)), ...
    min(min(plotcn(end),plotgm(end)),plotrb(end)), 1.01]);     
xlabel('Iterations')
ylabel('Relative Nonlinear Residuals')
legend('red-black','newt-gmres','newton')

  disp(' ')
  disp(' *** final results ***')
  disp(' ')
  disp('         |    red-black    |    newt-gmres    |    classic      |   ')
fprintf('      it |      %3i        |      %3i         |     %3i         | \n', ...
    itrb-1,itgm-1,it-1)
fprintf('    keff |%16.13f |%16.13f  |%16.13f | \n', x1(end),x2(end-1),x3(end-1))

n2=norm( feval('respfctAUG',x2) );
n3=norm( feval('respfctAUG',x3) );
fprintf('||F(x)|| |    %4.3e   |    %4.3e    |    %4.3e   | \n',n1,n2,n3)
fprintf('       p |      %4.3f      |      %4.3f       |      %4.3f      | \n',p1,p2,p3)
fprintf('    time |      %4.3f      |      %4.3f       |      %4.3f      | \n',t1,t2,t3)


%profile viewer